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O'Dea et al. (1983, J. Phys. Chem. 97, 391 1-3918) proposed an empirical procedure for obtaining estimates 
and confidence intervals for kinetic parameters in a model for pulse voltammetric data. Their goal was to find 
a procedure that would run in real time, not necessarily one that would have well-defined statistical properties. 
In this paper we investigate some of the statistical properties of their procedure. We show that their estimation 
method is equivalent to maximum likelihood estimation, and their confidence intervals, while related to like- 
lihood ratio confidence regions, have a coverage probability that is not fixed and that is potentially quite large. 
We suggest modifications of their procedure that lead to more traditional confidence intervals. We examine the 
effect on their procedure of the presence of nuisance paramters. Finally we discuss the possibility of serially 
correlated errors. 
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1. Introduction 

O'Dea et al. (1983) proposed a nonlinear regression 
procedure for estimating, and obtaining confidence in- 
tervals for, kinetic parameters describing the reduction 
of Zn(II) at a stationary mercury electrode in aqueous 
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solutions of NaN0 3 . In this paper we examine the statis- 
tical properties of the procedure and suggest mod- 
ifications to improve these properties. 

In section 2 we describe O'Dea's procedure. In sec- 
tion 3 we show his estimation procedure to be equiv- 
alent to maximum likelihood estimation. In section 4 we 
show his interval estimation procedure produces inter- 
vals that are related to higher-dimensional confidence 
regions obtained by likelihood ratio theory, and we sug- 
gest modifications to the procedure that will produce 
confidence intervals with the desired coverage proba- 
bility. In section 5 we question the assumption of 
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independent errors and examine the effect of including 
another parameter in the model to describe the apparent 
autoregressive error structure. 

The notation used here is similar to that used by 
O'Dea, but as is customary in literature on regression we 
use Y as the dependent variable. 



2. Description of the Procedure 

O'Dea models the observed response at time t\ to an 
arbitrary pulse sequence by 

Y t =afi+c+e t , 

where a and c are unknown constants that convey no 
kinetic information, {ej is a sequence of errors that are 
assumed to be independent with mean and unknown 
variance cr 2 . The function f\=f{t x , a, k ,Ej) is the solution 
of an integral equation. It depends on unknown kinetic 
parameters a, k, and Ej, and it must be obtained by 
solving the integral equation numerically. 

The kinetic parameters are of primary interest, so 
O'Dea uses a nonlinear optimization procedure to find 
the values of these parameters that maximize the cor- 
relation R between Y and f(t ,a,k ,Ej). These values are 
taken as the estimates. Estimates of a and c can then be 
obtained by simple linear regression of Y on f(t,a,k,Ej), 
and o- can be estimated by the standard deviation of the 
residuals from this regression. 

With no error the correlation R calculated above 
would be equal to unity. O'Dea measures the deviation 
from unity by R = 1 — R and defines R min as the optimum 
value of R . To measure the uncertainty in his estimate of 
a he fixes k and Ej at their optimal values and finds the 
two values of a that give R = 3 R mio . He calls the interval 
between these values a "confidence interval" for a, but 
he assigns no confidence level to the interval. He com- 
putes similar intervals for k and Ej. 



3. Maximum Likelihood Estimation 

A more traditional approach to a problem of this sort 
would be to write the likelihood or log likelihood func- 
tion for the problem and maximize it as a function of the 
unknown parameters. For normally distributed errors 
this is equivalent to choosing parameter values that min- 
imze the sum of squared residuals. In this section we 
show that the above estimation procedure is also equiv- 
alent to maximum likelihood estimation. 

If n is the number of observations, the log likelihood 
L is given by 



L(a,c, a,k,Ej,(T)= -- logger 2 )- 



■is 



Y-c-afi 



The maximum is easily found by noting that for any 
value of cr, the expression is maximized by choosing a, 
c, a, k, and Ej to minimize the sum of squared residuals 
SS R =1(Yi—c —afy 1 . It is a simple matter to show 
that if Y+n-^Yi and SS T = J,(Y i -Y) 1 , then 
SS R =(l—R 2 )SSr> so SSk is a minimum when R is a 
maximum (in absolute value). Therefore O'Dea's esti- 
mates are the maximum likelihood estimates. 

4. Confidence Intervals 

Confidence regions for unknown parameters are often 
found by computing the maximum likelihood estimates 
and then finding other sets of parameter values for 
which the likelihood function, or an approximation to 
the likelihood function, is not much smaller. O'Dea's 
procedure is related to this approach. 

Define L (a,k,Ej) as the maximum over a, c, and cr of 
the log likelihood L(a,c,a,k,Ej,cr). Using the re- 
lationships between R, SSr, and SS T above and max- 
imizing over or gives 

Z(o,^)=(-«/2)(l+Iog(2fl-) + logS5 r 

+ log(l-i? 2 )-log«). 

O'Dea's procedure involves finding six points — the two 
endpoints of the confidence intervals for each 
parameter — with the same correlation R = I — 3R min . 
The above expression for L(a,k,Ej) shows that these 
points have the same log likelihood as well. 

Let a,k~, and Ex be the maximum likelihood estimates 
of a, k, and Ej. The quantity X = ex.p(L(a,/e,Ej) 
—L(a,k,Ej)) is called the likelihood ratio. It can be 
shown that 2 log X has an asymptotic chi-square distri- 
bution with 3 degrees of freedom if a, k, and Ej are the 
true parameter values. Then P[2 log \<7.815]=0.95, so 
the parameter values for which 2 log \<7.815 form a 
95% confidence region for the true values of the param- 
eters. This region is bounded by the roughly ellipsoidal 
surface A(a,£,2?J)=exp(3.908), as shown in figure 1. In 
general, any surface of constant X bounds some con- 
fidence region. 

The confidence level of the region bounded by the 
surface containing O'Dea's points can be determined by 
computing X. In his procedure 

21og?t=-«[log(l-{l-^ ram ) 2 )~log(l-~(l-3^ mla ) 2 )] 
;s -nlog(2ii min /6 J R min )=Hlog3, 

since (^ n ,i n ) 2 «^ mi , 1 . Here re =81, so 2logk~ 89. By 
comparison, the region bounded by the surface for 
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Figure 1-Relationship between O'Dea's intervals and a three- 
dimensional confidence ellipsoid. The ellipsoid passes through the 
endpoints of the solid line segments. 

which 21ogX= 16.268 has a confidence level of 99.9%, 
so a three-dimensional confidence region found using 
21og\=89 would be very conservative. 

A more customary confidence level is 95%. Since the 
likelihood ratio can be written as a function of the cor- 
relation, it is possible to use a modification of O'Dea's 
procedure to find points on the boundary of a 95% 
confidence region. Rather than increasing R by a factor 
of 3, the appropriate factor is the value of b for which 81 
log b =7.815, or b zs 1.10. For example, in a sample data 
set that does not appear in O'Dea's original paper, 
a = .22522. Increasing R by a factor of 3 produces the 
interval (.22139, .22916], while the factor 1.10 leads to 
the interval [.22435, .22609]. 

On the other hand O'Dea's goal was not to find points 
in a confidence region for all three parameters, but to 
find separate confidence intervals for each parameter. In 
order to use the distribution of the likelihood ratio, 
O'Dea's procedure must be modified so that in com- 
puting the endpoints of the confidence interval for one 
parameter, the likelihood is maximized over the other 
two parameters. Twice the log of this likelihood ratio 
has an asymptotic chi-square distribution with one de- 
gree of freedom. 

In the case of a, for example, this is done by com- 
paring 21ogk=2[L(a,£,£i)-L(a,!c(cL),E±(ay\ to a chi- 
square distribution with one degree of freedom, where 
k~(s) and Z?i(s) are the values of k and E± that maximize 
L (p.,k,E^) subject to the restriction a=s. The 95% point 
of this distribution is 3.841, so the values of a for which 
2 log \<3.841 form a 95% confidence interval for the 
true parameter value. 

This is best illustrated in two dimensions, as in figure 
2. Here approximately elliptical contours of constant X 
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Figure 2-Comparison of confidence interval with interval computed 
by O'Dea's procedure in two dimensions. 

are plotted as a function of a and k for constant E±. (The 
complete contours are ellipsoids in three dimensions.) 
The inner ellipse has 21ogX= 3.841, while the outer el- 
lipse has 21ogX=7.815. The endpoints of the confidence 
interval for a are the points on the contour that have 
tangents perpendicular to the a axis. This interval can be 
compared with the interval found by O'Dea's pro- 
cedure, which is that portion of the k= films that is 
within the outer ellipse. 

Which is larger? If the two ellipses have major axes 
parallel to the coordinate axes, O'Dea's intervals are 
longer and their coverage probabilities exceed 95%. 
While this is not desirable, it increases the probability 
that his interval will contain the true parameters. But if 
the major axes are not parallel to the coordinate axes and 
if the lengths of the minor axes are small, O'Dea's inter- 
vals are shorter and have a coverage probability less 
than 95%. Unfortunately it is not possible to determine 
which is the case by looking only at the points examined 
in his procedure. 

There are two sensible remedies to this problem. The 
first, the likelihood ratio method, is similar in spirit to 
O'Dea's original procedure. This method involves find- 
ing the confidence interval as described above by find- 
ing those values of the first parameter that produce the 
proper likelihood ratio when the likelihood is max- 
imized over the other two parameters. In this problem, 
though, it is time consuming to calculate/; and its deriv- 
atives do not have simple analytic expressions, so re- 
peated maximization of the likelihood may be too com- 
putationally burdensome. 

The other method, an asymptotic normal approxi- 
mation, is the one we use here. This involves assuming 
that the maximum likelihood estimates have a multi- 
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variate normal distribution with a mean vector equal to 
the true parameter values and with covariance matrix 
equal to minus the inverse of the second derivative of 
the log likelihood L . (This is equivalent to the likelihood 
ratio method applied to a quadratic approximation to 
the log likelihood.) Since there is no analytic expression 
for the second derivative in this problem, we use a nu- 
merical approximation. 

In the example given above, the estimated covariance 
matrix is 



S = 



1.0334 


-.6908 


.0396 


.6908 


8.6984 


-.7509 


.0396 


-.7509 


.1384 



xio- 



A 95% confidence interval for a is given by 
a±1.96(5„)i, or [.22459, .22585]. This is narrower than 
the interval obtained above using O'Dea's procedure 
with a factor of 1.10. 



5. Residual Autocorrelation 

The above derivations are valid if the errors {ej are 
independent normal random variables with a common 
variance. In practice this assumption must be checked. 
This is especially true when, as in this case, mea- 
surements are taken over time. It is often reasonable to 



suspect that measurements at neighboring time points 
may be correlated. 

The errors {ej are not observed, but they can be 
estimated by the residuals, or the differences between 
the observed Y { and the fitted values 
?i=c+af(ti,a,h£x). Figure 3 is a plot of Y andf as a 
function of time. Figure 4 is a plot of the residuals 
e\ = Yi— Y-, over time. If the residuals were independent 
we would not expect to find any pattern here, but in fact 
there is a pronounced tendency for residuals at neigh- 
boring time points to have the same sign. 

There are three possible causes for this phenomenon. 
First, it is possible that this is an artifact of the fitting 
procedure. Even when the errors are independent, fit- 
ting an ordinary linear regression produces residuals 
that have some correlation (for example, they sum to 
zero). It is possible that minimizing the sum of squared 
residuals in this more complicated model produces re- 
siduals with some autocorrelation. However experi- 
ments with the model do not support this hypothesis. 

A second possible cause is model inadequacy. The 
model relates an imposed voltage to an observed cur- 
rent, and the voltage is highly correlated with time. If 
the equation used here is not the true relationship be- 
tween the current and the voltage, there may be cor- 
relation between the current and the residuals, and this 
dependence could be masquerading as time dependence. 
The cure for this difficulty is to propose alternative 
models that better fit the data. 




Figure 3-Observed 
Current. 



and Fitted 
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Figure 4-Residuals as a Function 
of Time. 
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The third possible cause is actual autocorrelation in 
the errors, and this autocorrelation can be modeled as 
well. We proceed under the assumption that the true 
errors have time-dependent correlation. 

Two common models for time series are the first or- 
der autoregressive model 

ei=pei_i-t-Wi 

and the first order moving average model 

e i =M i -|-pa i _i, 

where in both cases {«J is a sequence of independent 
normal random variables with mean and common un- 
known variance, and p is an unknown parameter be- 
tween — 1 and 1. Other possible models are the higher 
order models, where terms from earlier time points are 
used, and mixed models, where e, is modeled as a linear 
combination of ei-i,...,Gi- p and «j,...,Hj_ ? . 

Two tools useful for identifying a good model are the 
autocorrelation function and the partial autocorrelation 
function. These appear in figures 5 and 6. The sample 
autocorrelation function is simply the correlation of e, 
and e\_ k plotted as a function of k. For a moving average 
process of order q the true autocorrelation function is 
for k >q. For autoregressive processes and mixed pro- 
cesses the true autocorrelation function approaches as 
k— >oo, but it is not identically for all k beyond some 
finite value. The sample autocorrelation function in fig- 



ure 5 seems to be more consistent with that of the auto- 
regressive and mixed models, since there does not seem 
to be a sharp cutoff. 

The partial autocorrelation function is more compli- 
cated, but its interpretation is quite simple. It is the 
"dual" of the autocorrelation function, in that it is 
for all k >p for an autoregressive process of order/? , and 
it approaches as k—*oo but it does not vanish for mov- 
ing average and mixed processes. The sample partial 
autocorrelation function in figure 6 shows a large value 
at k — 1 and smaller values for k > 1. It is never exactly 
0, but for most k values the sample partial auto- 
correlation falls inside the boundary that marks the val- 
ues that are significantly different from 0. The function 
seems to be consistent with what might be expected 
from a first order autoregressive process. 

The new model 



with 



£i=pei_i+Hi 



is equivalent to 

r,=py 1 _,+a(/I-p/_i)+c(l-p)-l-«i, 

which is a nonlinear regression model with independent 
errors. 

There are now four parameters to be estimated, in 
addition to a, c, and cr. But if only the original three 
parameters are of interest, it is possible to treat p as one 
of the nuisance parameters by using a variant of the 
Cochrane-Orcutt procedure, as follows: 
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Figure 5-RestduaI Autocorrela- 
tion Function. 
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Figure 6-Residual Partial Auto- 
correlation Function. 
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1) For any given a, k, and E{, compute {ft}. 

2) Estimate a and c by linear regression to get {e}. 

3) Estimate p by the sample correlation of the {e}. 

4) Regress Y-pY,^ on f-pf\-i to get new esti- 
mates of a and c; and new residuals {e}. 

5) Repeat steps 3 and 4 until convergence. 

6) Compute the sum of squares 2«?=2(ej— pei_i) 2 . 

The computer time needed for these steps is much less 



than that needed to compute {/!}, so the estimation is 
much faster if the nonlinear optimization program 
searches only in the three-dimensional space of {a,k,E±), 
For each set of trial parameter values the above steps 
can be performed to minimize the residual sum of 
squares over the nuisance parameters. The resulting esti- 
mate of a is .22473. 

The other calculation can also be repeated for this 
new model. The estimated covariance matrix is 
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s = 



4.8344 


-5.6256 


.4766 


5.6256 


37.5982 


-2.6193 


.4766 


-2.6193 


.6513 



xio- 



terval that is about twice as long as the one in the inde- 
pendence model, but still only a third as long as the 
interval obtained by O'Dea's procedure. 



This is roughly four times the previous covariance ma- 
trix, so the length of the new confidence interval is 
about twice that of the previous confidence interval. 
The new interval is [.22336, .22600]. 

The four intervals around a are compared in figure 7. 
The procedure used in O'Dea's original paper produces 
an interval obtained from a three dimensional con- 
fidence ellipsoid with a very large confidence level, 
and it is quite long. The interval is shortened by using a 
95% confidence ellipsoid, but it still does not have a 
95% coverage probability. The 95% confidence inter- 
val is still shorter. Taking into account the apparent 
autoregressive error structure leads to a confidence in- 
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